This knitr document contains the code used to analyze and visualize data from the VIS camera system and water data. Knit-ing this document will generate an HTML document that contains both the embedded R code and the output.
Attach the required R packages
library(qvalue, warn.conflicts = FALSE)
library(ggplot2, warn.conflicts = FALSE)
library(lubridate, warn.conflicts = FALSE)
library(MASS, warn.conflicts = FALSE)
library(car, warn.conflicts = FALSE)
## Warning: package 'car' was built under R version 3.1.2
library(nlme, warn.conflicts = FALSE)
library(mvtnorm, warn.conflicts = FALSE)
library(grid, warn.conflicts = FALSE)
Use linear regression to create a simple model for using water volume to predict soil volume water content.
Download data for soil wet/dry weight and volume water content measurements.
if (!file.exists('soil_weigth_vwc.txt')) {
download.file('http://files.figshare.com/1939954/soil_weigth_vwc.txt',
'soil_weigth_vwc.txt')
}
Read the soil volume water content data
vwc.data = read.table(file="soil_weigth_vwc.txt", sep="\t", header=TRUE)
Calculate the average soil dry weight
mean(vwc.data$weight_dry)
## [1] 72.92138
Create a linear model for water volume and volume water content. Water volumes >= 260 appear to have saturated the soil water carrying capacity, so remove them from the model.
vwc.lm = lm(vwc_wet ~ water_vol, data=vwc.data[vwc.data$water_vol < 260,])
summary(vwc.lm)
##
## Call:
## lm(formula = vwc_wet ~ water_vol, data = vwc.data[vwc.data$water_vol <
## 260, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.5426 -1.6196 0.3976 1.4214 4.0771
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.032863 1.389638 -2.182 0.0391 *
## water_vol 0.233197 0.008145 28.630 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.082 on 24 degrees of freedom
## Multiple R-squared: 0.9716, Adjusted R-squared: 0.9704
## F-statistic: 819.7 on 1 and 24 DF, p-value: < 2.2e-16
Plot the model
vwc.plot = ggplot(vwc.data[vwc.data$water_vol < 260,], aes(x=water_vol, y=vwc_wet)) +
geom_point(size=2) +
geom_smooth(method='lm', formula=y~x) +
scale_x_continuous("Water volume (ml)") +
scale_y_continuous("Volume water content (%)") +
theme_bw() +
theme(axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold"))
print(vwc.plot)
Predict volume water contents for the water treatment groups
treatment.water = data.frame(water_vol=c(217, 144.5, 72))
treatment.water$vwc = predict(object = vwc.lm, newdata=treatment.water)
print(treatment.water)
## water_vol vwc
## 1 217.0 47.57084
## 2 144.5 30.66407
## 3 72.0 13.75731
LemnaTec VIS camera zoom units range from 1 to 6000, which correspond to 1 to 6X zoom.
zoom.lm = lm(zoom.camera ~ zoom, data=data.frame(zoom=c(1,6000),
zoom.camera=c(1,6)))
summary(zoom.lm)
##
## Call:
## lm(formula = zoom.camera ~ zoom, data = data.frame(zoom = c(1,
## 6000), zoom.camera = c(1, 6)))
##
## Residuals:
## ALL 2 residuals are 0: no residual degrees of freedom!
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.9991665 NA NA NA
## zoom 0.0008335 NA NA NA
##
## Residual standard error: NaN on 0 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: NaN
## F-statistic: NaN on 1 and 0 DF, p-value: NA
In this section we define models that are used to convert area and length between camera zoom levels to a common scale.
Download data for a reference object imaged at different zoom levels.
if (!file.exists('zoom_calibration_data.txt')) {
download.file('http://files.figshare.com/1845355/zoom_calibration_data.txt',
'zoom_calibration_data.txt')
}
Read zoom calibrartion data
z.data = read.table(file="zoom_calibration_data.txt", sep="\t", header=TRUE)
Calculate px per cm
z.data$px_cm = z.data$length_px / z.data$length_cm
Convert LemnaTec zoom units to camera zoom units
z.data$zoom.camera = predict(object = zoom.lm, newdata=z.data)
Fit a variety of regression models to relative object area by zoom level. Test whether an exponential or 2nd order polynomial model fits best (lowest AIC)
Non-linear regression (exponential)
area.coef = coef(nls(log(rel_area) ~ log(a * exp(b * zoom.camera)),
z.data, start = c(a = 1, b = 0.01)))
area.coef = data.frame(a=area.coef[1], b=area.coef[2])
area.nls = nls(rel_area ~ a * exp(b * zoom.camera),
data = z.data, start=c(a=area.coef$a, b=area.coef$b))
summary(area.nls)
##
## Formula: rel_area ~ a * exp(b * zoom.camera)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## a 0.37245 0.02159 17.25 <2e-16 ***
## b 0.86668 0.01504 57.61 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3306 on 32 degrees of freedom
##
## Number of iterations to convergence: 3
## Achieved convergence tolerance: 1.605e-06
Non-linear regression (polynomial)
area.pol = lm(rel_area ~ zoom.camera + I(zoom.camera^2), z.data)
summary(area.pol)
##
## Call:
## lm(formula = rel_area ~ zoom.camera + I(zoom.camera^2), data = z.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.8494 -0.3820 0.1176 0.3173 0.7494
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.42793 0.48203 9.186 2.34e-10 ***
## zoom.camera -4.68583 0.40917 -11.452 1.14e-12 ***
## I(zoom.camera^2) 1.64781 0.07786 21.163 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.431 on 31 degrees of freedom
## Multiple R-squared: 0.9907, Adjusted R-squared: 0.99
## F-statistic: 1643 on 2 and 31 DF, p-value: < 2.2e-16
AIC
AIC(area.nls, area.pol)
## df AIC
## area.nls 3 25.16270
## area.pol 4 44.12291
The exponential model minimizes AIC. Plot exponential model.
nls.plot = ggplot(z.data, aes(x=zoom.camera, y=rel_area)) +
geom_point(size=2.5) +
scale_x_continuous(lim=c(0.5,4.5), "Camera zoom setting") +
scale_y_continuous(lim=c(0,16),
"Reference object relative pixel area") +
stat_smooth(data=z.data, aes(x=zoom.camera, y=rel_area),
method="nls", se=FALSE, formula=y ~ a * exp(b * x),
start=c(a=area.coef$a, b=area.coef$b)) +
theme_bw() +
theme(axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold"))
nls.plot = nls.plot + labs(title='Figure 12A')
print(nls.plot)
Fit a variety of regression models to px/cm by zoom level. Test whether an exponential or 2nd order polynomial model fits best (lowest AIC).
Non-linear regression (exponential)
len.coef = coef(nls(log(px_cm) ~ log(a * exp(b * zoom.camera)),
z.data[z.data$camera == 'VIS SV',], start = c(a = 1, b = 0.01)))
len.coef = data.frame(a=len.coef[1], b=len.coef[2])
len.nls = nls(px_cm ~ a * exp(b * zoom.camera),
data = z.data[z.data$camera == 'VIS SV',],
start=c(a=len.coef$a, b=len.coef$b))
summary(len.nls)
##
## Formula: px_cm ~ a * exp(b * zoom.camera)
##
## Parameters:
## Estimate Std. Error t value Pr(>|t|)
## a 8.370121 0.150430 55.64 <2e-16 ***
## b 0.445769 0.005504 80.98 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5106 on 14 degrees of freedom
##
## Number of iterations to convergence: 3
## Achieved convergence tolerance: 2.15e-06
Length zoom correction using a 2 order polynomial (px/cm given a zoom setting). Length correction only works for side-view images right now because scale in top-view images are affected by the plant lifter position in addition to zoom.
len.poly = lm(px_cm ~ zoom.camera + I(zoom.camera^2),
data=z.data[z.data$camera == 'VIS SV',])
summary(len.poly)
##
## Call:
## lm(formula = px_cm ~ zoom.camera + I(zoom.camera^2), data = z.data[z.data$camera ==
## "VIS SV", ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.29376 -0.11153 0.02836 0.14223 0.21043
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.89815 0.30540 48.78 4.13e-16 ***
## zoom.camera -3.74980 0.27285 -13.74 4.04e-09 ***
## I(zoom.camera^2) 3.12322 0.05472 57.08 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1741 on 13 degrees of freedom
## Multiple R-squared: 0.9998, Adjusted R-squared: 0.9998
## F-statistic: 3.402e+04 on 2 and 13 DF, p-value: < 2.2e-16
AIC
AIC(len.nls, len.poly)
## df AIC
## len.nls 3 27.762784
## len.poly 4 -5.857059
The polynomial model minimizes AIC. Plot polynomial model.
poly.plot = ggplot(z.data[z.data$camera == 'VIS SV',], aes(x=zoom.camera, y=px_cm)) +
geom_point(size=2.5) +
scale_x_continuous(lim=c(0.5,4.5), "Camera zoom setting") +
scale_y_continuous(lim=c(0,50),
"Reference object length scale (px/cm)") +
stat_smooth(data=z.data[z.data$camera == 'VIS SV',],
aes(x=zoom.camera, y=px_cm), method="lm",
formula=y ~ x + I(x^2)) +
theme_bw() +
theme(axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold"))
poly.plot = poly.plot + labs(title='Figure 12B')
print(poly.plot)
In this section we measure how well PlantCV estimates plant height.
Download data manually measured plant height data set (n = 173).
if (!file.exists('height_model_data.txt')) {
download.file('http://files.figshare.com/1845361/height_model_data.txt',
'height_model_data.txt')
}
Read height data.
ht.data = read.table(file="height_model_data.txt",sep="\t",header=TRUE)
Convert LemnaTec zoom units to camera zoom units
ht.data$zoom.camera = predict(object = zoom.lm, newdata=ht.data)
Convert pixel traits to cm with zoom correction. px/cm calibration
px_cm = predict(object = len.poly, newdata=ht.data)
ht.data$height_above_bound_cm = ht.data$height_above_bound / px_cm
Height linear model
height.model = lm(manual_height~height_above_bound_cm, ht.data)
summary(height.model)
##
## Call:
## lm(formula = manual_height ~ height_above_bound_cm, data = ht.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.9229 -0.3034 -0.1263 0.3916 5.2015
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.303813 0.149776 2.028 0.0441 *
## height_above_bound_cm 0.852158 0.002789 305.528 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.103 on 171 degrees of freedom
## Multiple R-squared: 0.9982, Adjusted R-squared: 0.9982
## F-statistic: 9.335e+04 on 1 and 171 DF, p-value: < 2.2e-16
Plot the height linear model.
hlm.plot = ggplot(ht.data, aes(x=height_above_bound_cm, y=manual_height)) +
geom_point(aes(color=factor(round(zoom.camera,1))),size=2.5) +
geom_smooth(method="lm", formula=y~x, color='black') +
scale_x_continuous(lim=c(0,110), breaks=c(0,25,50,75,100), "Estimated height (cm)") +
scale_y_continuous(lim=c(0,110), breaks=c(0,25,50,75,100), "Manually measured height (cm)") +
theme_bw() +
theme(legend.position=c(0.2,0.75),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Zoom")
hlm.plot = hlm.plot + labs(title='Figure 3A')
print(hlm.plot)
In this section we model fresh- and dry-weight above ground biomass using image measurements.
Download data manually measured plant biomass data set (n = 41).
if (!file.exists('biomass_model_data.txt')) {
download.file('http://files.figshare.com/1845360/biomass_model_data.txt',
'biomass_model_data.txt')
}
Read biomass data.
st.data = read.table(file='biomass_model_data.txt', sep="\t", header=TRUE,
stringsAsFactors=FALSE)
Create out-of-frame indicator variable.
st.data$outind = NA
st.data[st.data$out_of_frame == T,]$outind = 1
st.data[st.data$out_of_frame == F,]$outind = 0
Genotype indicator variable.
st.data$group = NA
st.data[st.data$genotype == 'A10',]$group = 0
st.data[st.data$genotype == 'B100',]$group = 1
Full model. Includes side-view area, top-view area and height.
fw.full = lm(fresh_weight ~ sv_area * tv_area * height_above_bound, st.data)
Step-wise model selection with AIC.
fw.step = stepAIC(fw.full, direction="both")
## Start: AIC=-50.4
## fresh_weight ~ sv_area * tv_area * height_above_bound
##
## Df Sum of Sq RSS AIC
## - sv_area:tv_area:height_above_bound 1 0.055428 8.1731 -52.122
## <none> 8.1177 -50.401
##
## Step: AIC=-52.12
## fresh_weight ~ sv_area + tv_area + height_above_bound + sv_area:tv_area +
## sv_area:height_above_bound + tv_area:height_above_bound
##
## Df Sum of Sq RSS AIC
## - tv_area:height_above_bound 1 0.024335 8.1974 -54.000
## - sv_area:tv_area 1 0.061951 8.2350 -53.812
## - sv_area:height_above_bound 1 0.108786 8.2819 -53.580
## <none> 8.1731 -52.122
## + sv_area:tv_area:height_above_bound 1 0.055428 8.1177 -50.401
##
## Step: AIC=-54
## fresh_weight ~ sv_area + tv_area + height_above_bound + sv_area:tv_area +
## sv_area:height_above_bound
##
## Df Sum of Sq RSS AIC
## - sv_area:tv_area 1 0.31153 8.5089 -54.471
## <none> 8.1974 -54.000
## - sv_area:height_above_bound 1 0.68121 8.8786 -52.727
## + tv_area:height_above_bound 1 0.02433 8.1731 -52.122
##
## Step: AIC=-54.47
## fresh_weight ~ sv_area + tv_area + height_above_bound + sv_area:height_above_bound
##
## Df Sum of Sq RSS AIC
## <none> 8.5089 -54.471
## - sv_area:height_above_bound 1 0.42578 8.9347 -54.469
## + sv_area:tv_area 1 0.31153 8.1974 -54.000
## - tv_area 1 0.54733 9.0563 -53.915
## + tv_area:height_above_bound 1 0.27391 8.2350 -53.812
summary(fw.step)
##
## Call:
## lm(formula = fresh_weight ~ sv_area + tv_area + height_above_bound +
## sv_area:height_above_bound, data = st.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.19711 -0.10790 -0.00974 0.08153 1.55179
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.885e-01 2.513e-01 0.750 0.4581
## sv_area 4.673e-05 4.812e-06 9.710 1.36e-11 ***
## tv_area -1.045e-05 6.864e-06 -1.522 0.1368
## height_above_bound -3.280e-02 1.451e-02 -2.260 0.0299 *
## sv_area:height_above_bound 1.166e-07 8.688e-08 1.342 0.1879
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4862 on 36 degrees of freedom
## Multiple R-squared: 0.9835, Adjusted R-squared: 0.9816
## F-statistic: 535.9 on 4 and 36 DF, p-value: < 2.2e-16
AIC model
fw.aic = lm(fresh_weight ~ sv_area + tv_area + height_above_bound +
sv_area*height_above_bound, st.data)
summary(fw.aic)
##
## Call:
## lm(formula = fresh_weight ~ sv_area + tv_area + height_above_bound +
## sv_area * height_above_bound, data = st.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.19711 -0.10790 -0.00974 0.08153 1.55179
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.885e-01 2.513e-01 0.750 0.4581
## sv_area 4.673e-05 4.812e-06 9.710 1.36e-11 ***
## tv_area -1.045e-05 6.864e-06 -1.522 0.1368
## height_above_bound -3.280e-02 1.451e-02 -2.260 0.0299 *
## sv_area:height_above_bound 1.166e-07 8.688e-08 1.342 0.1879
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4862 on 36 degrees of freedom
## Multiple R-squared: 0.9835, Adjusted R-squared: 0.9816
## F-statistic: 535.9 on 4 and 36 DF, p-value: < 2.2e-16
The AIC model contains tv_area and height which does not have a significant coefficient, test dropping.
fw.red = lm(fresh_weight ~ sv_area, st.data)
summary(fw.red)
##
## Call:
## lm(formula = fresh_weight ~ sv_area, data = st.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.17949 -0.25179 0.02035 0.23804 1.28235
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.044e-01 1.060e-01 -2.872 0.00657 **
## sv_area 3.935e-05 8.757e-07 44.931 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5003 on 39 degrees of freedom
## Multiple R-squared: 0.981, Adjusted R-squared: 0.9806
## F-statistic: 2019 on 1 and 39 DF, p-value: < 2.2e-16
Goodness of fit.
anova(fw.aic, fw.red)
## Analysis of Variance Table
##
## Model 1: fresh_weight ~ sv_area + tv_area + height_above_bound + sv_area *
## height_above_bound
## Model 2: fresh_weight ~ sv_area
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 36 8.5089
## 2 39 9.7630 -3 -1.2541 1.7686 0.1706
SV area model.
sv.model = lm(fresh_weight ~ sv_area, st.data)
summary(sv.model)
##
## Call:
## lm(formula = fresh_weight ~ sv_area, data = st.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.17949 -0.25179 0.02035 0.23804 1.28235
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.044e-01 1.060e-01 -2.872 0.00657 **
## sv_area 3.935e-05 8.757e-07 44.931 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5003 on 39 degrees of freedom
## Multiple R-squared: 0.981, Adjusted R-squared: 0.9806
## F-statistic: 2019 on 1 and 39 DF, p-value: < 2.2e-16
Plot SV model
sv.model.plot = ggplot(st.data,aes(x=sv_area/1e5, y=fresh_weight)) +
geom_smooth(method="lm", color="black", formula = y ~ x) +
geom_point(size=2.5) +
scale_x_continuous("Shoot and leaf area (x10^5 px)") +
scale_y_continuous("Fresh-weight biomass (g)") +
theme_bw() +
theme(axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold"))
sv.model.plot = sv.model.plot + labs(title='Figure 4A')
print(sv.model.plot)
SV area model with out-of-frame.
sv.ind.model = lm(fresh_weight ~ sv_area + outind, st.data)
anova(sv.model, sv.ind.model)
## Analysis of Variance Table
##
## Model 1: fresh_weight ~ sv_area
## Model 2: fresh_weight ~ sv_area + outind
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 39 9.7630
## 2 38 7.5498 1 2.2132 11.139 0.0019 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(sv.ind.model)
##
## Call:
## lm(formula = fresh_weight ~ sv_area + outind, data = st.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.75671 -0.21287 -0.00372 0.22304 1.17010
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.499e-01 9.583e-02 -2.607 0.0130 *
## sv_area 4.028e-05 8.288e-07 48.599 <2e-16 ***
## outind -5.963e-01 1.787e-01 -3.338 0.0019 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4457 on 38 degrees of freedom
## Multiple R-squared: 0.9853, Adjusted R-squared: 0.9846
## F-statistic: 1277 on 2 and 38 DF, p-value: < 2.2e-16
Plot SV model with out-of-frame.
sv.model.out.plot = ggplot(st.data, aes(x=sv_area/1e5, y=fresh_weight,
group=out_of_frame, color=out_of_frame)) +
geom_point(size=2.5) +
geom_smooth(method="lm", color="black") +
scale_x_continuous("Shoot and leaf area (x10^5 px)") +
scale_y_continuous("Fresh-weight biomass (g)") +
theme_bw() +
theme(legend.position=c(0.2,0.8),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Out-of-frame")
print(sv.model.out.plot)
SV area model with genotype
sv.gt.model = lm(fresh_weight ~ sv_area + group, st.data)
summary(sv.gt.model)
##
## Call:
## lm(formula = fresh_weight ~ sv_area + group, data = st.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.32109 -0.29602 0.03093 0.21960 1.15107
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.510e-01 1.342e-01 -3.362 0.00178 **
## sv_area 3.956e-05 8.637e-07 45.802 < 2e-16 ***
## group 2.647e-01 1.542e-01 1.717 0.09419 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4883 on 38 degrees of freedom
## Multiple R-squared: 0.9824, Adjusted R-squared: 0.9815
## F-statistic: 1061 on 2 and 38 DF, p-value: < 2.2e-16
Plot SV model with genotype.
sv.model.gen.plot = ggplot(st.data, aes(x=sv_area/1e5, y=fresh_weight,
group=genotype, color=genotype)) +
geom_point(size=2.5) +
geom_smooth(method="lm", color="black") +
scale_x_continuous("Shoot and leaf area (x10^5 px)") +
scale_y_continuous("Fresh-weight biomass (g)") +
theme_bw() +
theme(legend.position=c(0.2,0.8),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Genotype")
print(sv.model.gen.plot)
SV area model
dry.sv.model = lm(dry_weight ~ sv_area, st.data)
summary(dry.sv.model)
##
## Call:
## lm(formula = dry_weight ~ sv_area, data = st.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.183938 -0.047916 0.009925 0.039315 0.195700
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.295e-02 1.479e-02 -3.58 0.000938 ***
## sv_area 4.634e-06 1.222e-07 37.93 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06981 on 39 degrees of freedom
## Multiple R-squared: 0.9736, Adjusted R-squared: 0.9729
## F-statistic: 1438 on 1 and 39 DF, p-value: < 2.2e-16
Plot dry-weight side-view model
dry.sv.model.plot = ggplot(st.data,aes(x=sv_area/1e5, y=dry_weight)) +
geom_smooth(method="lm", color="black", formula = y ~ x) +
geom_point(size=2.5) +
scale_x_continuous("Shoot and leaf area (x10^5 px)") +
scale_y_continuous("Dry-weight biomass (g)") +
theme_bw() +
theme(axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold"))
dry.sv.model.plot = dry.sv.model.plot + labs(title='Supplemental Figure S2')
print(dry.sv.model.plot)
SV area model with out-of-frame.
dry.sv.ind.model = lm(dry_weight ~ sv_area + outind, st.data)
summary(dry.sv.ind.model)
##
## Call:
## lm(formula = dry_weight ~ sv_area + outind, data = st.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.133525 -0.027886 0.006673 0.024790 0.131393
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.393e-02 1.255e-02 -3.501 0.001202 **
## sv_area 4.789e-06 1.085e-07 44.124 < 2e-16 ***
## outind -9.869e-02 2.340e-02 -4.219 0.000147 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05836 on 38 degrees of freedom
## Multiple R-squared: 0.982, Adjusted R-squared: 0.9811
## F-statistic: 1038 on 2 and 38 DF, p-value: < 2.2e-16
dry.sv.model.out.plot = ggplot(st.data, aes(x=sv_area/1e5, y=dry_weight,
group=out_of_frame, color=out_of_frame)) +
geom_point(size=2.5) +
geom_smooth(method="lm", color="black") +
scale_x_continuous("Shoot and leaf area (x10^5 px)") +
scale_y_continuous("Dry-weight biomass (g)") +
theme_bw() +
theme(legend.position=c(0.2,0.8),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Out-of-frame")
print(dry.sv.model.out.plot)
SV area model with genotypes
dry.sv.gt.model = lm(dry_weight ~ sv_area + group, st.data)
summary(dry.sv.gt.model)
##
## Call:
## lm(formula = dry_weight ~ sv_area + group, data = st.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.183197 -0.049205 0.008654 0.040345 0.196263
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.413e-02 1.943e-02 -2.786 0.00828 **
## sv_area 4.636e-06 1.251e-07 37.061 < 2e-16 ***
## group 2.130e-03 2.233e-02 0.095 0.92451
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07071 on 38 degrees of freedom
## Multiple R-squared: 0.9736, Adjusted R-squared: 0.9722
## F-statistic: 701 on 2 and 38 DF, p-value: < 2.2e-16
dry.sv.model.gen.plot = ggplot(st.data, aes(x=sv_area/1e5, y=dry_weight,
group=genotype, color=genotype)) +
geom_point(size=2.5) +
geom_smooth(method="lm", color="black") +
scale_x_continuous("Shoot and leaf area (x10^5 px)") +
scale_y_continuous("Dry-weight biomass (g)") +
theme_bw() +
theme(legend.position=c(0.2,0.8),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Genotype")
print(dry.sv.model.gen.plot)
In this section we analyze geometric traits from the complete VIS data set.
Download the complete VIS data set. There were a total of 6,399 snapshots with VIS image data, but the download only includes the 6,207 snapshots that were successfully processed by PlantCV. Failed snapshots generally are those that lack a plant (empty pot controls, dead plants, etc.)
if (!file.exists('vis_snapshots.txt')) {
download.file('http://files.figshare.com/1845362/vis_snapshots.txt',
'vis_snapshots.txt')
}
Read data and format for analysis
vis.data = read.table(file="vis_snapshots.txt", sep=",", header=TRUE)
Planting date
planting_date = as.POSIXct("2013-11-26")
Add water treatment column coded in barcodes.
vis.data$treatment <- NA
vis.data$treatment[grep("AA", vis.data$plant_id)] <- 100
vis.data$treatment[grep("AB", vis.data$plant_id)] <- 0
vis.data$treatment[grep("AC", vis.data$plant_id)] <- 16
vis.data$treatment[grep("AD", vis.data$plant_id)] <- 33
vis.data$treatment[grep("AE", vis.data$plant_id)] <- 66
Add plant genotype column coded in barcodes.
vis.data$genotype <- NA
vis.data$genotype[grep("p1", vis.data$plant_id)] <- 'A10'
vis.data$genotype[grep("p2", vis.data$plant_id)] <- 'B100'
vis.data$genotype[grep("r1", vis.data$plant_id)] <- 'R20'
vis.data$genotype[grep("r2", vis.data$plant_id)] <- 'R70'
vis.data$genotype[grep("r3", vis.data$plant_id)] <- 'R98'
vis.data$genotype[grep("r4", vis.data$plant_id)] <- 'R102'
vis.data$genotype[grep("r5", vis.data$plant_id)] <- 'R128'
vis.data$genotype[grep("r6", vis.data$plant_id)] <- 'R133'
vis.data$genotype[grep("r7", vis.data$plant_id)] <- 'R161'
vis.data$genotype[grep("r8", vis.data$plant_id)] <- 'R187'
Add genotype x treatment group column
vis.data$group = paste(vis.data$genotype,'-',vis.data$treatment,sep='')
Remove plants that were sampled for biomass measurements.
vis.data = vis.data[!vis.data$plant_id %in% st.data$barcode,]
Add calendar-time data column using the Unix-time data
vis.data$date = as.POSIXct(vis.data$datetime, origin = "1970-01-01")
Calculate days after planting from planting data
vis.data$dap = as.numeric(vis.data$date - planting_date)
Predict fresh- and dry-weight biomass from linear models
vis.data$fw_biomass = predict.lm(object = sv.model, newdata=vis.data)
vis.data$dw_biomass = predict.lm(object = dry.sv.model, newdata=vis.data)
Plot height for S. viridis and S. italica water treatments 100% and 33% full-capacity.
height.plot = ggplot(vis.data[(vis.data$genotype == 'A10' |
vis.data$genotype == 'B100') &
(vis.data$treatment == 100 |
vis.data$treatment == 33),],
aes(x=dap, y=height_above_bound, color=factor(group))) +
geom_point(size=2.5) +
geom_smooth(method="loess",size=1) +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(lim=c(0,120),
breaks=c(0,20,40,60,80,100,120),
name="Estimated height (cm)") +
theme_bw() +
theme(legend.position=c(0.25,0.75),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Genotype-treatment")
height.plot = height.plot + labs(title='Figure 3B')
print(height.plot)
Plot height for all genotypes.
height.facet.plot = ggplot(vis.data[vis.data$treatment == 100 |
vis.data$treatment == 33,],
aes(x=dap, y=height_above_bound, color=factor(treatment))) +
geom_point(size=1) +
geom_smooth(method="loess",size=1) +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(lim=c(0,120),
breaks=c(0,20,40,60,80,100,120),
name="Estimated height (cm)") +
theme_bw() +
theme(legend.position=c(0.8,0.1),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Treatment") +
facet_wrap(~ genotype, ncol = 3)
height.facet.plot = height.facet.plot + labs(title='Supplemental Figure S1')
print(height.facet.plot)
For each day, calculate the mean and 95% confidence intervals for each genotype in the control water group.
height_per_day = function(vis.data) {
dap = c()
genotype = c()
int.low = c()
int.up = c()
est = c()
# Loop through each day and genotype
for(d in min(as.integer(vis.data$dap)):max(as.integer(vis.data$dap))) {
if(d != 24) {
for(g in levels(factor(vis.data$genotype))) {
h.data = vis.data[vis.data$genotype == g & vis.data$treatment == 100 &
as.integer(vis.data$dap) == d,]$height_above_bound
if(sd(h.data) != 0) {
dap = c(dap,d)
genotype = c(genotype,g)
test = t.test(h.data)
int.low = c(int.low,test$conf.int[1])
int.up = c(int.up,test$conf.int[2])
est = c(est,test$estimate)
}
}
}
}
#drought_resp = drought_resp[-14,]
results = data.frame(dap=as.numeric(dap),
genotype=genotype,
conf.int.low=int.low,
conf.int.up=int.up,
mean=est)
return(results)
}
Height per day 95% CI. Write results to file.
height.perday.results = height_per_day(vis.data)
write.csv(height.perday.results,file="height.perday.100water.csv")
print(height.perday.results)
## dap genotype conf.int.low conf.int.up mean
## 1 11 A10 7.8384656 9.927418 8.882942
## 2 11 B100 5.3761414 6.677522 6.026831
## 3 11 R102 2.6078989 6.747923 4.677911
## 4 11 R128 3.4255116 6.098742 4.762127
## 5 11 R133 5.8106341 11.190037 8.500336
## 6 11 R161 5.9550014 6.742792 6.348897
## 7 11 R187 4.6001185 8.238236 6.419177
## 8 11 R20 8.7685119 12.487009 10.627761
## 9 11 R70 6.9589765 8.860251 7.909614
## 10 11 R98 11.3958077 13.499331 12.447569
## 11 12 A10 10.5706419 14.040344 12.305493
## 12 12 B100 7.1779957 9.543168 8.360582
## 13 12 R102 4.1386602 10.394107 7.266383
## 14 12 R128 4.7644938 9.358302 7.061398
## 15 12 R133 0.3769178 32.979175 16.678046
## 16 12 R161 5.5667113 10.395653 7.981182
## 17 12 R187 8.9117000 14.337900 11.624800
## 18 12 R20 -24.8325684 57.616115 16.391774
## 19 12 R70 6.2768196 11.896208 9.086514
## 20 12 R98 14.4614239 19.097635 16.779529
## 21 13 A10 14.2175454 16.845714 15.531630
## 22 13 B100 10.9868413 12.771688 11.879265
## 23 13 R102 8.1006537 12.638227 10.369440
## 24 13 R128 10.0067593 13.083564 11.545162
## 25 13 R133 15.6105089 20.434429 18.022469
## 26 13 R161 11.5326626 14.465157 12.998910
## 27 13 R187 11.0472672 15.088387 13.067827
## 28 13 R20 14.8393242 21.736455 18.287890
## 29 13 R70 12.7527374 14.743597 13.748167
## 30 13 R98 17.9034635 24.295281 21.099372
## 31 14 A10 16.0353867 20.806689 18.421038
## 32 14 B100 14.8758324 16.482712 15.679272
## 33 14 R102 9.5848380 17.897359 13.741099
## 34 14 R128 11.5441914 18.807803 15.175997
## 35 14 R133 24.4394863 26.317757 25.378622
## 36 14 R161 11.3351757 16.721924 14.028550
## 37 14 R187 15.9765377 22.156071 19.066305
## 38 14 R20 9.1279616 35.869895 22.498929
## 39 14 R70 12.8683843 21.043290 16.955837
## 40 14 R98 22.7445324 26.529747 24.637140
## 41 15 A10 20.7800983 23.377751 22.078924
## 42 15 B100 17.7777752 20.518164 19.147970
## 43 15 R102 11.0573174 17.922030 14.489674
## 44 15 R128 16.7501844 21.295836 19.023010
## 45 15 R133 27.0393189 33.948830 30.494074
## 46 15 R161 17.1755797 20.612348 18.893964
## 47 15 R187 15.8387658 22.954208 19.396487
## 48 15 R20 25.1788370 33.253535 29.216186
## 49 15 R70 18.9601854 22.863978 20.912081
## 50 15 R98 29.0498999 31.955600 30.502750
## 51 16 A10 20.0577561 25.806257 22.932007
## 52 16 B100 20.0413990 23.203604 21.622501
## 53 16 R102 14.8670912 21.144546 18.005818
## 54 16 R128 16.5582461 30.333804 23.446025
## 55 16 R133 33.0362095 35.768486 34.402348
## 56 16 R161 16.4543384 23.558591 20.006465
## 57 16 R187 22.3958994 28.061383 25.228641
## 58 16 R20 9.0056082 50.912573 29.959091
## 59 16 R70 20.2903256 31.959948 26.125137
## 60 16 R98 32.7722175 35.627625 34.199921
## 61 17 A10 27.4881214 32.177027 29.832574
## 62 17 B100 25.7375286 29.145299 27.441414
## 63 17 R102 9.3832005 41.668709 25.525955
## 64 17 R128 26.8195976 35.438632 31.129115
## 65 17 R133 35.0134141 35.046324 35.029869
## 66 17 R161 18.5118930 41.274711 29.893302
## 67 17 R187 26.7251684 31.553360 29.139264
## 68 17 R20 34.5958844 35.180457 34.888171
## 69 17 R70 27.8333370 32.092941 29.963139
## 70 17 R98 34.9761994 35.063296 35.019748
## 71 18 A10 25.4649891 34.113739 29.789364
## 72 18 B100 29.5995162 33.060621 31.330069
## 73 18 R102 23.9458634 35.780013 29.862938
## 74 18 R128 28.8943167 33.278893 31.086605
## 75 18 R161 22.7950766 33.274113 28.034595
## 76 18 R187 28.6642353 38.039275 33.351755
## 77 18 R20 34.8319033 35.217714 35.024809
## 78 18 R70 32.5635770 36.332210 34.447894
## 79 18 R98 35.0377193 35.072626 35.055172
## 80 19 A10 29.9565044 33.715470 31.835987
## 81 19 B100 28.7587854 32.088308 30.423547
## 82 19 R102 18.6004422 36.868409 27.734426
## 83 19 R128 35.0349633 35.062368 35.048666
## 84 19 R133 35.0369049 35.055222 35.046063
## 85 19 R161 27.5528428 35.523189 31.538016
## 86 19 R187 30.2853896 36.052235 33.168812
## 87 19 R20 35.0431394 35.064445 35.053792
## 88 19 R70 32.8273407 35.931809 34.379575
## 89 20 A10 30.7762330 35.328096 33.052164
## 90 20 B100 31.4095265 35.224147 33.316837
## 91 20 R102 29.8675676 36.244858 33.056213
## 92 20 R128 31.8182665 36.308301 34.063284
## 93 20 R133 35.0292025 35.062924 35.046063
## 94 20 R161 27.6114760 34.513152 31.062314
## 95 20 R187 33.7013663 35.619516 34.660441
## 96 20 R20 35.0065633 35.093660 35.050112
## 97 20 R70 31.1659064 36.873618 34.019762
## 98 20 R98 35.0336568 35.066567 35.050112
## 99 21 A10 34.7350643 43.619676 39.177370
## 100 21 B100 33.9185376 40.511408 37.214973
## 101 21 R102 20.4335385 53.700974 37.067256
## 102 21 R128 45.7878729 59.113583 52.450728
## 103 21 R133 69.5748771 75.552885 72.563881
## 104 21 R161 30.8383653 46.826065 38.832215
## 105 21 R187 41.2621776 55.259613 48.260895
## 106 21 R20 56.3428640 68.941224 62.642044
## 107 21 R70 45.2148076 52.242602 48.728705
## 108 21 R98 56.5053136 69.874481 63.189898
## 109 22 A10 42.1333456 51.279316 46.706331
## 110 22 B100 35.1894676 46.590960 40.890214
## 111 22 R102 32.4270279 54.660958 43.543993
## 112 22 R128 42.7607416 57.866538 50.313640
## 113 22 R133 64.7341342 95.014087 79.874111
## 114 22 R161 31.9080496 39.897530 35.902790
## 115 22 R187 49.2083816 68.125731 58.667056
## 116 22 R20 45.5047607 113.876490 79.690625
## 117 22 R70 37.0178551 70.222627 53.620241
## 118 22 R98 64.6229237 77.916294 71.269609
## 119 23 A10 46.1205983 57.372009 51.746303
## 120 23 B100 39.0566070 49.213411 44.135009
## 121 23 R102 22.1355965 62.113193 42.124395
## 122 23 R128 56.6189532 68.402234 62.510593
## 123 23 R133 80.1898844 87.496497 83.843190
## 124 23 R161 39.5976979 48.909886 44.253792
## 125 23 R187 44.9050613 74.178923 59.541992
## 126 23 R20 65.9802051 85.598085 75.789145
## 127 23 R70 58.6470448 67.651642 63.149343
## 128 23 R98 60.6980969 79.556520 70.127309
## 129 25 A10 55.9173941 66.531485 61.224439
## 130 25 B100 42.1436759 60.516430 51.330053
## 131 25 R102 37.5833391 74.825638 56.204488
## 132 25 R128 59.0014536 77.608288 68.304871
## 133 25 R133 91.5153533 106.509863 99.012608
## 134 25 R161 41.7599259 51.527936 46.643931
## 135 25 R187 69.4178315 80.320776 74.869304
## 136 25 R20 47.3113464 137.062553 92.186950
## 137 25 R70 54.8053892 93.180449 73.992919
## 138 25 R98 79.0174138 95.332386 87.174900
## 139 26 A10 60.7309054 74.484346 67.607626
## 140 26 B100 51.9573893 67.031311 59.494350
## 141 26 R102 29.0283221 90.704621 59.866472
## 142 26 R128 69.7432021 85.949025 77.846114
## 143 26 R133 102.1131687 102.209599 102.161384
## 144 26 R161 44.7659963 77.522223 61.144110
## 145 26 R187 64.2493428 83.643787 73.946565
## 146 26 R20 95.4818290 103.551591 99.516710
## 147 26 R70 80.8950037 90.615600 85.755302
## 148 26 R98 67.3596660 96.850149 82.104907
## 149 27 A10 68.1673027 79.966027 74.066665
## 150 27 B100 52.3930899 73.531038 62.962064
## 151 27 R102 46.4307081 91.299260 68.864984
## 152 27 R128 69.7163940 92.330224 81.023309
## 153 27 R133 102.0534520 102.623594 102.338523
## 154 27 R161 47.0038080 59.572263 53.288036
## 155 27 R187 79.9731196 87.133834 83.553477
## 156 27 R20 94.6611877 107.250057 100.955622
## 157 27 R70 67.7411474 102.378830 85.059989
## 158 27 R98 92.2825030 102.888066 97.585284
## 159 28 A10 72.4019654 84.816088 78.609027
## 160 28 B100 60.2668834 76.201220 68.234052
## 161 28 R102 38.3458163 100.670481 69.508149
## 162 28 R128 88.7635776 97.020262 92.891920
## 163 28 R133 102.0427075 102.412501 102.227604
## 164 28 R161 54.3065643 86.175685 70.241125
## 165 28 R187 68.2614296 91.127546 79.694488
## 166 28 R20 101.9335479 102.697144 102.315346
## 167 28 R70 96.9145446 101.216624 99.065584
## 168 28 R98 81.5682665 106.133203 93.850735
## 169 29 A10 77.9931551 87.733707 82.863431
## 170 29 B100 64.8732783 90.167093 77.520186
## 171 29 R102 51.4886408 103.817301 77.652971
## 172 29 R128 78.2635038 103.522285 90.892894
## 173 29 R133 102.1483036 102.853222 102.500763
## 174 29 R161 60.8247415 72.598147 66.711444
## 175 29 R187 85.3621398 92.131055 88.746598
## 176 29 R20 102.1549445 102.730696 102.442820
## 177 29 R70 81.5101446 109.747356 95.628750
## 178 29 R98 100.9210632 102.477379 101.699221
## 179 30 A10 79.7975021 93.165534 86.481518
## 180 30 B100 69.9817638 86.831826 78.406795
## 181 30 R102 52.6392582 105.567685 79.103472
## 182 30 R128 100.2996586 102.191014 101.245336
## 183 30 R133 102.1362130 102.517656 102.326935
## 184 30 R161 72.2477821 89.354608 80.801195
## 185 30 R187 75.6166966 95.754845 85.685771
## 186 30 R20 102.3143461 102.675591 102.494969
## 187 30 R70 102.0282615 102.642163 102.335212
## 188 30 R98 93.7598548 105.298398 99.529126
## 189 31 A10 87.1296524 94.989746 91.059699
## 190 31 B100 76.9270063 98.740991 87.833999
## 191 31 R102 62.8067625 104.995504 83.901133
## 192 31 R128 89.1628356 105.911162 97.536999
## 193 31 R133 102.3840680 102.756520 102.570294
## 194 31 R161 68.6018081 80.383543 74.492676
## 195 31 R187 89.8228715 96.405196 93.114034
## 196 31 R20 101.8550011 103.378296 102.616648
## 197 31 R70 87.6302155 109.560625 98.595420
## 198 31 R98 102.2831228 102.641146 102.462134
## 199 32 A10 86.6320427 101.920801 94.276422
## 200 32 B100 82.8496624 97.181019 90.015341
## 201 32 R102 66.9055230 110.225409 88.565466
## 202 32 R128 105.0325661 107.488859 106.260713
## 203 32 R133 107.1520311 107.382340 107.267186
## 204 32 R161 77.4568438 96.766731 87.111787
## 205 32 R187 85.8563435 97.724316 91.790330
## 206 32 R20 107.2913801 107.379954 107.335667
## 207 32 R70 107.0598209 107.278888 107.169355
## 208 32 R98 105.1210435 108.008120 106.564582
## 209 33 A10 96.0947018 103.735921 99.915311
## 210 33 B100 89.2731829 105.354814 97.313999
## 211 33 R102 75.8742679 110.271342 93.072805
## 212 33 R128 96.1577692 110.345571 103.251670
## 213 33 R133 107.1725115 107.361860 107.267186
## 214 33 R161 81.8858866 88.986237 85.436062
## 215 33 R187 93.8716625 99.386701 96.629182
## 216 33 R20 107.0615755 107.597308 107.329442
## 217 33 R70 92.3612921 114.279011 103.320152
## 218 33 R98 107.1300106 107.383609 107.256810
Statistical analysis for height differences. Use pairwise comparisons on S. viridis and S. italica over two-day intervals.
analyze_height = function(vis.data, genotype) {
days = c()
diff.low = c()
diff.up = c()
pvals = c()
for(day in levels(factor(as.integer(vis.data$dap)))) {
day = as.integer(day)
control = vis.data[(as.integer(vis.data$dap) == day |
as.integer(vis.data$dap) == day + 1) &
vis.data$genotype == genotype &
vis.data$treatment == 100,]$height_above_bound
drought = vis.data[(as.integer(vis.data$dap) == day |
as.integer(vis.data$dap) == day + 1) &
vis.data$genotype == genotype &
vis.data$treatment == 33,]$height_above_bound
test = t.test(x=control, y=drought)
days = c(days, day)
diff.low = c(diff.low, test$conf.int[1])
diff.up = c(diff.up, test$conf.int[2])
pvals = c(pvals, test$p.value)
}
results = data.frame(dap=as.numeric(days),
conf.int.low=diff.low,
conf.int.up=diff.up,
pvalue=pvals)
return(results)
}
Test for treatment effect on two-day intervals
a10.height.results = analyze_height(vis.data, 'A10')
b100.height.results = analyze_height(vis.data, 'B100')
Control for multiple testing by controlling the FDR
qvalues.height = qvalue(c(a10.height.results$pvalue, b100.height.results$pvalue))
a10.height.results$qvalue = qvalues.height$qvalues[1:nrow(a10.height.results)]
b100.height.results$qvalue = qvalues.height$qvalues[
(nrow(a10.height.results)+1):(nrow(a10.height.results) + nrow(b100.height.results))]
Assign genotypes for merged table
a10.height.results$genotype = 'A10'
b100.height.results$genotype = 'B100'
print(a10.height.results)
## dap conf.int.low conf.int.up pvalue qvalue genotype
## 1 11 1.3757087 5.139721 1.622327e-03 1.057500e-03 A10
## 2 12 -0.7380484 4.515406 1.490756e-01 5.492423e-02 A10
## 3 13 -0.5912495 3.543282 1.562446e-01 5.516696e-02 A10
## 4 14 -1.1319577 4.329657 2.381965e-01 6.960212e-02 A10
## 5 15 -0.1299653 5.167933 6.143652e-02 3.062407e-02 A10
## 6 16 -1.3607720 7.080051 1.715524e-01 5.519138e-02 A10
## 7 17 -2.9248819 5.400151 5.353411e-01 1.191457e-01 A10
## 8 18 -0.4959663 6.673738 8.769306e-02 3.911081e-02 A10
## 9 19 1.4409925 7.969806 6.845578e-03 4.143499e-03 A10
## 10 20 3.6149867 12.029236 6.790664e-04 4.795304e-04 A10
## 11 21 6.7803624 17.104831 5.343934e-05 6.694424e-05 A10
## 12 22 10.2910019 21.102986 1.471517e-06 1.246954e-05 A10
## 13 23 9.2657053 24.761796 2.155551e-04 1.660545e-04 A10
## 14 25 7.2790852 20.484417 1.699507e-04 1.600168e-04 A10
## 15 26 8.5545706 23.844353 1.911306e-04 1.619628e-04 A10
## 16 27 11.5264141 26.176502 1.946187e-05 3.758465e-05 A10
## 17 28 10.7457158 25.976630 5.606520e-05 6.694424e-05 A10
## 18 29 11.4473342 24.918441 7.632409e-06 3.233827e-05 A10
## 19 30 11.9787099 27.114832 2.217662e-05 3.758465e-05 A10
## 20 31 12.9567838 28.606319 1.247030e-05 3.522418e-05 A10
## 21 32 12.2483513 29.889421 6.320013e-05 6.694424e-05 A10
## 22 33 3.4267550 38.584764 2.690657e-02 1.518525e-02 A10
print(b100.height.results)
## dap conf.int.low conf.int.up pvalue qvalue genotype
## 1 11 -2.1888257 0.4486175 0.18236609 0.05519138 B100
## 2 12 -2.2019131 1.9489490 0.89939162 0.17724157 B100
## 3 13 -3.7808622 0.5052064 0.12627973 0.05095650 B100
## 4 14 -2.7496601 2.5699624 0.94357005 0.18172165 B100
## 5 15 -4.8854099 0.8999123 0.16277335 0.05517324 B100
## 6 16 -2.8974999 3.9598806 0.74490405 0.15029213 B100
## 7 17 -5.0026064 2.9536026 0.58132070 0.12014816 B100
## 8 18 -1.9644558 5.1388196 0.35534030 0.09713326 B100
## 9 19 -2.3699180 4.9749886 0.46324842 0.10904271 B100
## 10 20 -0.8187179 6.7951985 0.11836288 0.05014997 B100
## 11 21 -8.9817009 4.6032442 0.50430459 0.11549851 B100
## 12 22 -8.3643015 4.7370840 0.57202885 0.12014816 B100
## 13 23 -11.0131272 5.4448269 0.46303807 0.10904271 B100
## 14 25 -13.4045213 2.5892801 0.17774764 0.05519138 B100
## 15 26 -17.3140056 -1.0219101 0.02867192 0.01518525 B100
## 16 27 -14.3161656 2.1808693 0.14381616 0.05492423 B100
## 17 28 -17.0011012 1.2317531 0.08754828 0.03911081 B100
## 18 29 -11.9653330 5.2377910 0.43167803 0.10904271 B100
## 19 30 -12.3111302 4.5263257 0.35381611 0.09713326 B100
## 20 31 -4.5854432 10.1681899 0.44686134 0.10904271 B100
## 21 32 -5.5221448 10.1865601 0.54834969 0.11914565 B100
## 22 33 -6.1404245 14.8562724 0.38830762 0.10282794 B100
Output the A10 and B100 tables to the same file
write.table(a10.height.results, file='height.stats.csv', sep = ',',
row.names = FALSE, append = FALSE)
write.table(b100.height.results, file='height.stats.csv', sep = ',',
row.names = FALSE, append = TRUE, col.names = FALSE)
Plot fresh-weight biomass for S. viridis and S. italica water treatments 100% and 33% full-capacity.
biomass.plot = ggplot(vis.data[(vis.data$genotype == 'A10' |
vis.data$genotype == 'B100') &
(vis.data$treatment == 100 |
vis.data$treatment == 33),],
aes(x=dap, y=fw_biomass, color=factor(group))) +
geom_point(size=2.5) +
geom_smooth(method="loess",size=1) +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(lim=c(-1,41),
name="Estimated biomass (g)") +
theme_bw() +
theme(legend.position=c(0.25,0.75),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Genotype-treatment")
biomass.plot = biomass.plot + labs(title='Figure 4B')
print(biomass.plot)
Plot fresh-weight biomass for all genotypes.
biomass.facet.plot = ggplot(vis.data[vis.data$treatment == 100 |
vis.data$treatment == 33,],
aes(x=dap, y=fw_biomass, color=factor(treatment))) +
geom_point(size=1) +
geom_smooth(method="loess",size=1) +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(lim=c(-1,41),
name="Estimated biomass (g)") +
theme_bw() +
theme(legend.position=c(0.8,0.1),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Treatment") +
facet_wrap(~ genotype, ncol = 3)
biomass.facet.plot = biomass.facet.plot + labs(title='Supplemental Figure S3')
print(biomass.facet.plot)
Statistical analysis for fresh-weight biomass differences. Use pairwise comparisons on S. viridis and S. italica over two-day intervals.
analyze_biomass = function(vis.data, genotype) {
days = c()
diff.low = c()
diff.up = c()
pvals = c()
for(day in levels(factor(as.integer(vis.data$dap)))) {
day = as.integer(day)
control = vis.data[(as.integer(vis.data$dap) == day |
as.integer(vis.data$dap) == day + 1) &
vis.data$genotype == genotype &
vis.data$treatment == 100,]$fw_biomass
drought = vis.data[(as.integer(vis.data$dap) == day |
as.integer(vis.data$dap) == day + 1) &
vis.data$genotype == genotype &
vis.data$treatment == 33,]$fw_biomass
# If the control group has a lot more replicates (e.g. A10 and B100), randomly sample to drought sample size
if (genotype == 'A10' | genotype == 'B100') {
control = sample(control, size = length(drought))
}
test = t.test(x=control, y=drought)
days = c(days, day)
diff.low = c(diff.low, test$conf.int[1])
diff.up = c(diff.up, test$conf.int[2])
pvals = c(pvals, test$p.value)
}
results = data.frame(dap=as.numeric(days),
conf.int.low=diff.low,
conf.int.up=diff.up,
pvalue=pvals)
return(results)
}
Test for treatment effect on two-day intervals
a10.biomass.results = analyze_biomass(vis.data, 'A10')
b100.biomass.results = analyze_biomass(vis.data, 'B100')
Control for multiple testing by controlling the FDR
qvalues.biomass = qvalue(c(a10.biomass.results$pvalue, b100.biomass.results$pvalue))
a10.biomass.results$qvalue = qvalues.biomass$qvalues[1:nrow(a10.biomass.results)]
b100.biomass.results$qvalue = qvalues.biomass$qvalues[
(nrow(a10.biomass.results)+1):(nrow(a10.biomass.results) +
nrow(b100.biomass.results))]
Assign genotypes for merged table
a10.biomass.results$genotype = 'A10'
b100.biomass.results$genotype = 'B100'
print(a10.biomass.results)
## dap conf.int.low conf.int.up pvalue qvalue genotype
## 1 11 0.02914367 0.2201204 1.323444e-02 1.515363e-02 A10
## 2 12 -0.14788575 0.1756949 8.610677e-01 5.049914e-01 A10
## 3 13 -0.16165563 0.2625467 6.306207e-01 3.790869e-01 A10
## 4 14 -0.18405590 0.4646706 3.796492e-01 2.552294e-01 A10
## 5 15 -0.27354779 0.6885945 3.821226e-01 2.552294e-01 A10
## 6 16 -0.39537050 0.6661931 5.981374e-01 3.784843e-01 A10
## 7 17 0.10518290 3.1083996 3.780833e-02 3.952667e-02 A10
## 8 18 -0.13454475 2.2779082 7.835338e-02 7.138478e-02 A10
## 9 19 0.55639946 3.2418022 7.520995e-03 9.518142e-03 A10
## 10 20 1.79793954 4.9902898 2.723816e-04 4.093438e-04 A10
## 11 21 2.01188934 6.2302349 6.230113e-04 8.812062e-04 A10
## 12 22 4.48502027 8.5942426 2.003136e-06 4.013837e-06 A10
## 13 23 4.38451816 10.3131422 2.161645e-04 3.465163e-04 A10
## 14 25 6.34140762 10.9450081 1.001322e-07 3.009637e-07 A10
## 15 26 6.34064811 11.1993657 2.033517e-07 5.432952e-07 A10
## 16 27 8.43389988 12.1214602 4.849913e-11 2.332354e-10 A10
## 17 28 9.48931032 12.8236353 6.467846e-13 7.776071e-12 A10
## 18 29 8.43000027 11.5285952 1.588593e-12 1.041987e-11 A10
## 19 30 10.37389366 13.1633241 1.396552e-14 3.358054e-13 A10
## 20 31 9.37549132 12.8536239 1.733372e-12 1.041987e-11 A10
## 21 32 8.77790792 12.5759125 8.718795e-11 3.494103e-10 A10
## 22 33 9.13675701 14.3597517 1.794587e-06 4.013837e-06 A10
print(b100.biomass.results)
## dap conf.int.low conf.int.up pvalue qvalue genotype
## 1 11 -0.1298704 -0.02014007 9.958668e-03 1.197297e-02 B100
## 2 12 -0.0805700 0.14440699 5.635119e-01 3.662114e-01 B100
## 3 13 -0.2655002 0.05797019 1.994123e-01 1.453010e-01 B100
## 4 14 -0.2982824 0.27327839 9.285308e-01 5.315909e-01 B100
## 5 15 -0.5579089 0.03474439 8.015652e-02 7.138478e-02 B100
## 6 16 -0.4525040 0.74118118 6.175631e-01 3.790869e-01 B100
## 7 17 -1.0828481 1.08346063 9.995283e-01 5.462267e-01 B100
## 8 18 -0.2717842 1.41231582 1.737604e-01 1.347782e-01 B100
## 9 19 -1.1177500 1.15611705 9.725059e-01 5.438189e-01 B100
## 10 20 -0.2312011 2.65901500 9.508123e-02 8.165209e-02 B100
## 11 21 -0.2643127 2.63411000 1.034254e-01 8.575509e-02 B100
## 12 22 -0.5396748 3.20099463 1.535486e-01 1.230708e-01 B100
## 13 23 -2.9434030 7.26769789 3.233459e-01 2.286751e-01 B100
## 14 25 -1.2656478 5.61185394 1.982686e-01 1.453010e-01 B100
## 15 26 -0.1065401 5.98201981 5.765873e-02 5.545690e-02 B100
## 16 27 0.1611022 6.85503580 4.110506e-02 4.118268e-02 B100
## 17 28 0.7754796 9.21814913 2.356266e-02 2.575326e-02 B100
## 18 29 4.0254297 10.07486231 1.391536e-04 2.389995e-04 B100
## 19 30 6.9552625 11.09939185 7.806236e-09 2.681477e-08 B100
## 20 31 5.5617792 10.73510877 2.262272e-06 4.184388e-06 B100
## 21 32 3.0137106 11.18548970 2.080233e-03 2.778882e-03 B100
## 22 33 8.1298505 13.84599687 1.865143e-06 4.013837e-06 B100
Output the A10 and B100 tables to the same file
write.table(a10.biomass.results, file='biomass.stats.csv', sep = ',',
row.names = FALSE, append = FALSE)
write.table(b100.biomass.results, file='biomass.stats.csv', sep = ',',
row.names = FALSE, append = TRUE, col.names = FALSE)
Calculate the difference in biomass per genotype per day between the 33% and 100% water treatment groups.
drought_response = function(vis.data, genotypes) {
# Initialize drought response data frame with days
drought_resp = data.frame(day=c(
min(as.integer(vis.data$dap)):max(as.integer(vis.data$dap))))
# Initialize genotypes
for(g in genotypes) {
control = paste(g,'control',sep='')
drought = paste(g,'drought',sep='')
# Calculate median biomass per treatment per day
drought_resp[,paste(control)] = 0
drought_resp[,paste(drought)] = 0
for(d in min(as.integer(vis.data$dap)):max(as.integer(vis.data$dap))) {
drought_resp[drought_resp$day == d, paste(control)] =
median(vis.data[vis.data$genotype == g &
vis.data$treatment == 100 &
as.integer(vis.data$dap) == d,]$fw_biomass)
drought_resp[drought_resp$day == d, paste(drought)] =
median(vis.data[vis.data$genotype == g &
vis.data$treatment == 33 &
as.integer(vis.data$dap) == d,]$fw_biomass)
}
}
drought_resp = drought_resp[-14,]
# Calculate biomass loss to drought
days = c()
response = c()
genotype=c()
for(r in 1:nrow(drought_resp)) {
days = c(days, rep(drought_resp[r,]$day,length(genotypes)))
for(g in genotypes) {
control = paste(g,'control',sep='')
drought = paste(g,'drought',sep='')
response = c(response, drought_resp[r,paste(drought)] -
drought_resp[r,paste(control)])
genotype = c(genotype, g)
}
}
dr.df = data.frame(day=days,response=response,genotype=as.factor(genotype))
return(dr.df)
}
A10 vs B100 drought responses
drought.response.parents = drought_response(vis.data,c('A10','B100'))
Plot drought response results
genotype.colors = c('#E6A024','#4AB859')
resp.plot = ggplot(drought.response.parents,
aes(x=day, y=response, group=genotype, color=genotype)) +
geom_point(size=2.5) +
geom_smooth(method="loess") +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(name="Reduced accumulated biomass (g)") +
theme_bw() +
theme(legend.position=c(0.8,0.8),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
scale_colour_manual(values=genotype.colors) +
labs(color='Genotype')
resp.plot = resp.plot + labs(title='Figure 4C')
print(resp.plot)
In this experiment plants appear to follow relatively basic asymptotic growth patterns. Here we use non-linear least squares regression analysis to model a three-component logistic growth function for each genotype x treatment group.
First load three functions for non-linear growth curve analysis reproduced from:
Paine CET, Marthews TR, Vogt DR, Purves D, Rees M, Hector A, Turnbull LA (2012) How to fit nonlinear plant growth models and calculate growth rates: an update for ecologists. Methods in Ecology and Evolution 3: 245–256. doi: 10.1111/j.2041-210X.2011.00155.x.
output.logis.nlsList <- function(fit, times, CI = F, LOG = F, alpha = 0.05){
coef <- coef(fit)
params <- transform_param.logis(coef)
rates <- list()
groups <- rownames(params)
n.groups <- nrow(coef)
# compute rates for each group seperately
rates <- list()
for(i in 1:(n.groups)){
K <- params[i,1]; r <- params[i,2]; M0 <- params[i,3]
rates[[i]] = data.frame(
times = times,
M = (M0*K)/(M0+(K-M0)*exp(-r*times)),
AGR = (r*M0*K*(K-M0)*exp(-r*times))/(M0+(K-M0)*exp(-r*times))^2
)
rates[[i]]$RGRt <- rates[[i]]$AGR/rates[[i]]$M
rates[[i]]$RGRm <- r*(1 - rates[[i]]$M/K)
if(LOG == T){
rates[[i]]$RGRt <- rates[[i]]$AGR
rates[[i]]$RGRm <- r*rates[[i]]$M*(1-rates[[i]]$M/K)
rates[[i]]$AGR <- rates[[i]]$AGR*exp(rates[[i]]$M)
}
# commute CIs for each group's estaimates, if desired
if(CI == T){
cov <- summary(fit)$cov[i,,]
x <- y <- data.frame(rmvnorm(n=1000, mean=c(coef[i, 1], coef[i, 2], coef[i, 3]), sigma=cov))
x$K <- y[,1]
x$r <- 1/y[,3]
x$M0 <- y[,1]/(1 + exp(y[,2]/y[,3]))
M <- AGR <- RGRt <- RGRm <- matrix(NA, ncol = length(times), nrow = nrow(x))
for(j in 1:nrow(x)){
K <- x[j,4]; r <- x[j,5]; M0 <- x[j,6]
M[j,] <- (M0*K)/(M0+(K-M0)*exp(-r*times))
AGR[j,] <- (r*M0*K*(K-M0)*exp(-r*times))/(M0+(K-M0)*exp(-r*times))^2
RGRt[j,] <- AGR[j,]/M[j,]
RGRm[j,] <- r*(1 - M[j,]/K)
if(LOG ==T){
RGRt[j,] <- AGR[j,]
RGRm[j,] <- r*M[j,]*(1 - M[j,]/K)
AGR[j,] <- AGR[j,]*exp(M[j,])
}
}
}
CIs <- summarizer(list(M = M, AGR = AGR, RGRt = RGRt, RGRm = RGRm), alpha)
rates[[i]] <- cbind(rates[[i]], CIs)
}
names(rates) <- rownames(params)
# now compute differences among groups
diffs <- list()
for(i in 1:(n.groups-1)){
Ki <- params[i,1]; ri <- params[i,2]; M0i <- params[i,3]
for(j in (i+1):n.groups){
Kj <- params[j,1]; rj <- params[j,2]; M0j <- params[j,3]
diffs.ij = data.frame(
times = times,
diffM = rates[[i]]$M - rates[[j]]$M,
diffAGR = rates[[i]]$AGR - rates[[j]]$AGR,
diffRGRt = rates[[i]]$RGRt - rates[[j]]$RGRt
)
# comparing RGRm has to be done on a biomass basis. So it needs special treatment. First, we aev to know what range of biomasses are shared between groups.
Mmin <- max(min(rates[[i]]$M), min(rates[[j]]$M)) # yieds the range of overlapping masses between groups i & j
Mmax <- min(max(rates[[i]]$M), max(rates[[j]]$M))
diffs.ij$Mseq <- seq(Mmin, Mmax, length = 100)
if(LOG == F){
diffs.ij$diffRGRm <- ri*(1 - diffs.ij$Mseq/Ki) - rj*(1 - diffs.ij$Mseq/Kj)
} else{
diffs.ij$diffRGRm <- ri*diffs.ij$Mseq*(1 - diffs.ij$Mseq/Ki) - rj*Mseq*(1 - diffs.ij$Mseq/Kj)
}
}
if(CI == T){
# get params for group i
covi <- summary(fit)$cov[i,,]
xi <- yi <- data.frame(rmvnorm(n=1000, mean=c(coef[i, 1], coef[i, 2], coef[i, 3]), sigma=covi))
xi$K <- yi[,1]
xi$r <- 1/yi[,3]
xi$M0 <- yi[,1]/(1 + exp(yi[,2]/yi[,3]))
# get params for group j
covj <- summary(fit)$cov[j,,]
xj <- yj <- data.frame(rmvnorm(n=1000, mean=c(coef[j, 1], coef[j, 2], coef[j, 3]), sigma=covj))
xj$K <- yj[,1]
xj$r <- 1/yj[,3]
xj$M0 <- yj[,1]/(1 + exp(yj[,2]/yj[,3]))
# now compute diffs for each random set of drawn parameters
Mi <- Mj <- AGRi <- AGRj <- RGRti <- RGRtj <- RGRmi <- RGRmj <- diffM <- diffAGR <- diffRGRt <- diffRGRm <- matrix(NA, ncol = length(times), nrow = nrow(xi))
for(k in 1:nrow(xi)){
Ki <- xi[k,4]; ri <- xi[k,5]; M0i <- xi[k,6]
Kj <- xj[k,4]; rj <- xj[k,5]; M0j <- xj[k,6]
Mi[k,] <- (M0i*Ki)/(M0i+(Ki-M0i)*exp(-ri*times))
Mj[k,] <- (M0j*Kj)/(M0j+(Kj-M0j)*exp(-rj*times))
AGRi[k,] <- (ri*M0i*Ki*(Ki-M0i)*exp(-ri*times))/(M0i+(Ki-M0i)*exp(-ri*times))^2
AGRj[k,] <- (rj*M0j*Kj*(Kj-M0j)*exp(-rj*times))/(M0j+(Kj-M0j)*exp(-rj*times))^2
RGRti[k,] <- AGRi[k,]/Mi[k,]
RGRtj[k,] <- AGRj[k,]/Mj[k,]
RGRmi[k,] <- ri*(1 - diffs.ij$Mseq/Ki)
RGRmj[k,] <- rj*(1 - diffs.ij$Mseq/Kj)
if(LOG == T){
RGRti[k,] <- AGRi[k,]
RGRtj[k,] <- AGRj[k,]
RGRmi[k,] <- ri*diffs.ij$Mseq*(1 - diffs.ij$Mseq/Ki)
RGRmj[k,] <- rj*diffs.ij$Mseq*(1 - diffs.ij$Mseq/Kj)
AGRi[k,] <- AGRi[k,]*exp(Mi[k,])
AGRj[k,] <- AGRj[k,]*exp(Mj[k,])
}
diffM[k,] <- Mi[k,] - Mj[k,]
diffAGR[k,] <- AGRi[k,] - AGRj[k,]
diffRGRt[k,] <- RGRti[k,] - RGRtj[k,]
diffRGRm[k,] <- RGRmi[k,] - RGRmj[k,]
}
CIs <- summarizer(list(diffM = diffM, diffAGR = diffAGR, diffRGRt = diffRGRt, diffRGRm = diffRGRm), alpha)
diffs[[paste(groups[i], groups[j], sep = "_")]] <- cbind(diffs.ij, CIs)
} else{
diffs[[paste(groups[i], groups[j], sep = "_")]] <- diffs.ij
}
} # end loop over pairwise combinations of groups
out <- list(params = params, rates = rates, diffs = diffs)
return(out)
}
transform_param.logis <- function(coef){
K = coef[1]
r = 1/(coef[3])
M0 = K/(1 + exp(coef[2]/coef[3])) #untransform best-fit parameters to K, r and M0
if(is.data.frame(K)){
out <- cbind(K, r, M0)
} else {
out <- c(K, r, M0)
}
names(out) <- c("K", "r", "M0")
return(out)
}
# this function returns confidence envelopes around growth trajectories, and growth rates.
summarizer <- function(dat, alpha){
n <- length(dat)
quantiles <- c(alpha/2, 1-(alpha/2))
CIs <- data.frame(matrix(NA, ncol(dat[[1]]), n*2))
names(CIs) <- paste(rep(names(dat), each = 2), c("lo", "hi"), sep = ".")
for(i in 1:n){
CIs[,(2*i-1):(2*i)] <- t(apply(dat[[i]], 2, quantile, quantiles, na.rm = T))
}
return(CIs)
}
Subset the VIS data
sub.vis.data = vis.data[(vis.data$treatment == 100 | vis.data$treatment == 33),
c("plant_id","dap","fw_biomass","group")]
sub.vis.data$group = factor(sub.vis.data$group)
Group data
grouped.sub.vis.data = groupedData(fw_biomass ~ dap | group, sub.vis.data)
Fit three-component logistic functions for each genotype x treatment group
fit.logis = nlsList(fw_biomass ~ SSlogis(dap, Asym, xmid, scal),
data = grouped.sub.vis.data)
Output estimated growth rate parameters at even time intervals
est.interval = seq(min(sub.vis.data$dap), max(sub.vis.data$dap), length = 100)
out.fit.logis = output.logis.nlsList(fit.logis, times = est.interval,
CI = TRUE, LOG = FALSE, alpha = 0.05)
Compute the time and magnitude of maximum growth rate
logis.results = data.frame(group=levels(sub.vis.data$group))
logis.results$max.AGR = 0
logis.results$max.AGR.dap = 0
group.order = names(out.fit.logis$rates)
for(i in 1:length(group.order)) {
logis.results[logis.results == group.order[i],]$max.AGR.dap =
out.fit.logis$rates[[i]]$times[out.fit.logis$rates[[i]]$AGR ==
max(out.fit.logis$rates[[i]]$AGR)]
logis.results[logis.results == group.order[i],]$max.AGR =
max(out.fit.logis$rates[[i]]$AGR)
}
print(logis.results)
## group max.AGR max.AGR.dap
## 1 A10-100 2.994229 24.23556
## 2 A10-33 1.737080 25.36054
## 3 B100-100 2.389559 26.26053
## 4 B100-33 1.631987 24.68555
## 5 R102-100 2.593812 28.06050
## 6 R102-33 1.826986 26.48553
## 7 R128-100 2.995409 25.36054
## 8 R128-33 1.781943 24.23556
## 9 R133-100 3.267807 23.11057
## 10 R133-33 2.207510 23.33557
## 11 R161-100 2.858742 25.58554
## 12 R161-33 1.750648 24.68555
## 13 R187-100 2.444363 24.46055
## 14 R187-33 1.617136 24.46055
## 15 R20-100 3.240594 22.88557
## 16 R20-33 1.869087 22.66058
## 17 R70-100 2.878034 23.78556
## 18 R70-33 1.592505 23.33557
## 19 R98-100 2.509740 24.23556
## 20 R98-33 1.481989 24.68555
Plot logistic growth models for S. viridis and S. italica and absolute growth rates over time
parent.groups = c("A10-100", "A10-33", "B100-100", "B100-33")
group.colors = c("#F8766D", "#7CAE00", "#00BFC4", "#C77CFF")
# Biomass over time
gca.plot = ggplot(vis.data[(vis.data$genotype == 'A10' | vis.data$genotype == 'B100') &
(vis.data$treatment == 100 | vis.data$treatment == 33),],
aes(x=dap, y=fw_biomass, color=factor(group))) +
geom_point(size=2.5) +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(lim=c(-1,41), name="Estimated biomass (g)") +
theme_bw() +
theme(legend.position=c(0.2,0.8),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
scale_colour_manual(name='Genotype-treatment',
values=group.colors, labels=parent.groups)
agr.plot = ggplot() +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(name="Absolute growth rate (g/day)") +
theme_bw() +
theme(legend.position=c(0.2,0.8),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
scale_colour_manual(name='Genotype-treatment',
values=group.colors, labels=parent.groups)
# Add logistic growth models and confidence intervals
for(g in 1:length(parent.groups)) {
group = parent.groups[g]
i = match(group, group.order)
group.rates = as.data.frame(out.fit.logis$rates[i])
col.name = gsub('-', '.', group)
conf.int = data.frame(x=c(group.rates[,paste(col.name, '.times', sep='')],
rev(group.rates[,paste(col.name, '.times', sep='')])),
y=c(group.rates[,paste(col.name, '.M.lo', sep='')],
rev(group.rates[,paste(col.name, '.M.hi', sep='')])))
agr.conf.int = data.frame(x=c(group.rates[,paste(col.name, '.times', sep='')],
rev(group.rates[,paste(col.name, '.times', sep='')])),
y=c(group.rates[,paste(col.name, '.AGR.lo', sep='')],
rev(group.rates[,paste(col.name, '.AGR.hi', sep='')])))
gca.plot = gca.plot + geom_polygon(data=conf.int, aes(x=x, y=y),
fill='gray60', color=NA, alpha=0.4) +
geom_line(data=group.rates, aes_string(x=paste(col.name, '.times', sep=''),
y=paste(col.name, '.M', sep='')),
color=group.colors[g])
agr.plot = agr.plot +
geom_polygon(data=agr.conf.int, aes(x=x, y=y),
fill='gray60', color=NA, alpha=0.4) +
geom_line(data=group.rates, aes_string(x=paste(col.name, '.times', sep=''),
y=paste(col.name, '.AGR', sep='')),
color=group.colors[g])
}
gca.plot = gca.plot + labs(title='Figure 4B')
print(gca.plot)
agr.plot = agr.plot + labs(title='Figure 4C')
print(agr.plot)
Plot logistic growth models for all Setaria genotypes
plot_gca = function(genotype) {
groups = c(paste(genotype,'-100',sep=''), paste(genotype,'-33',sep=''))
group.colors = c("#00BFC4", "#F8766D")
# Biomass over time
gca.plot = ggplot(vis.data[(vis.data$genotype == genotype) &
(vis.data$treatment == 100 | vis.data$treatment == 33),],
aes(x=dap, y=fw_biomass, color=factor(treatment))) +
geom_point(size=1) +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(lim=c(-1,41), name="Estimated biomass (g)") +
theme_bw() +
theme(legend.position='none',
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Treatment") +
facet_grid(. ~ genotype)
# Add logistic growth models and confidence intervals
for(g in 1:length(groups)) {
group = groups[g]
i = match(group, group.order)
group.rates = as.data.frame(out.fit.logis$rates[i])
col.name = gsub('-', '.', group)
conf.int = data.frame(x=c(group.rates[,paste(col.name, '.times', sep='')],
rev(group.rates[,paste(col.name, '.times', sep='')])),
y=c(group.rates[,paste(col.name, '.M.lo', sep='')],
rev(group.rates[,paste(col.name, '.M.hi', sep='')])))
gca.plot = gca.plot + geom_polygon(data=conf.int, aes(x=x, y=y),
fill='gray60', color=NA, alpha=0.4) +
geom_line(data=group.rates, aes_string(x=paste(col.name, '.times', sep=''),
y=paste(col.name, '.M', sep='')),
color=group.colors[g])
}
return(gca.plot)
}
Multiple plot function modified from the Cookbook for R by Winston Chang: http://www.cookbook-r.com/Graphs/Multiple_graphs_on_one_page_(ggplot2)/
multiplot <- function(..., plotlist=NULL, cols=1, layout=NULL) {
# Make a list from the ... arguments and plotlist
plots <- c(list(...), plotlist)
numPlots = length(plots)
# If layout is NULL, then use 'cols' to determine layout
if (is.null(layout)) {
# Make the panel
# ncol: Number of columns of plots
# nrow: Number of rows needed, calculated from # of cols
layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
ncol = cols, nrow = ceiling(numPlots/cols),byrow=TRUE)
}
if (numPlots==1) {
print(plots[[1]])
} else {
# Set up the page
grid.newpage()
pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
# Make each plot, in the correct location
for (i in 1:numPlots) {
# Get the i,j matrix positions of the regions that contain this subplot
matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
layout.pos.col = matchidx$col))
}
}
}
Generate the Setaria biomass/growth plots
a10.gca.plot = plot_gca('A10')
b100.gca.plot = plot_gca('B100')
r102.gca.plot = plot_gca('R102')
r128.gca.plot = plot_gca('R128')
r133.gca.plot = plot_gca('R133')
r161.gca.plot = plot_gca('R161')
r187.gca.plot = plot_gca('R187')
r20.gca.plot = plot_gca('R20')
r70.gca.plot = plot_gca('R70')
r98.gca.plot = plot_gca('R98')
Supplemental Figure S3
multiplot(a10.gca.plot, b100.gca.plot, r102.gca.plot, r128.gca.plot,
r133.gca.plot, r161.gca.plot, r187.gca.plot, r20.gca.plot,
r70.gca.plot, r98.gca.plot, cols=3)
agr.plot = agr.plot + labs(title='Figure 4C')
print(agr.plot)
In this section we combine data on the volume of water added to each plant per water application (up to the target weight) with fresh-weight biomass data to estimate water use efficiency.
Download water volume data (n = 33,496 water applications).
if (!file.exists('B2_watering_data_phenofront.csv')) {
download.file('http://files.figshare.com/1845363/B2_watering_data_phenofront.csv',
'B2_watering_data_phenofront.csv')
}
Read data and format for analysis
water.data = read.table(file="B2_watering_data_phenofront.csv", sep=",", header=TRUE)
Remove the empty pot controls (Barcodes start with zero).
water.data = water.data[grep('000A', water.data$plant.barcode, invert = TRUE),]
Remove snapshots that were not valid waterings (water.amount = -1).
water.data = water.data[water.data$water.amount != -1,]
Add a treatment group column coded in the barcode.
water.data$treatment <- NA
water.data$treatment[grep("AA", water.data$plant.barcode)] <- 100
water.data$treatment[grep("AB", water.data$plant.barcode)] <- 0
water.data$treatment[grep("AC", water.data$plant.barcode)] <- 16
water.data$treatment[grep("AD", water.data$plant.barcode)] <- 33
water.data$treatment[grep("AE", water.data$plant.barcode)] <- 66
Add a plant genotype column coded in the barcode.
water.data$genotype <- NA
water.data$genotype[grep("p1", water.data$plant.barcode)] <- 'A10'
water.data$genotype[grep("p2", water.data$plant.barcode)] <- 'B100'
water.data$genotype[grep("r1", water.data$plant.barcode)] <- 'R20'
water.data$genotype[grep("r2", water.data$plant.barcode)] <- 'R70'
water.data$genotype[grep("r3", water.data$plant.barcode)] <- 'R98'
water.data$genotype[grep("r4", water.data$plant.barcode)] <- 'R102'
water.data$genotype[grep("r5", water.data$plant.barcode)] <- 'R128'
water.data$genotype[grep("r6", water.data$plant.barcode)] <- 'R133'
water.data$genotype[grep("r7", water.data$plant.barcode)] <- 'R161'
water.data$genotype[grep("r8", water.data$plant.barcode)] <- 'R187'
Add a genotype x treatment group column.
water.data$group = paste(water.data$genotype,'-',water.data$treatment,sep='')
Encode the calendar time column as a date-time.
water.data$timestamp = ymd_hms(water.data$timestamp)
Add a column for days after planting since the planting date.
water.data$dap = as.numeric(water.data$timestamp - planting_date)
Some plants were removed from the system but were not inactivated, so remove them.
bad.cars = unique(water.data$car.tag[water.data$water.amount >120 &
water.data$dap > 14])
water.data = water.data[!water.data$car.tag %in% bad.cars,]
Extract water data for S. viridis for Figure 1B.
sv.water = water.data[(water.data$treatment == 100 | water.data$treatment == 33) &
water.data$genotype == 'A10',]
Classify water job as early or later in the day. Early water treatments started at ZT23, later treatments started at ZT8.
sv.water$early.late <- NA #restrict to the scheduled watering later in the run
sv.water$early.late[hour(sv.water$timestamp) < 7] <- 'ZT23'
sv.water$early.late[hour(sv.water$timestamp) < 15 &
hour(sv.water$timestamp) > 11 ] <- 'ZT8'
sv.water = sv.water[!is.na(sv.water$early.late),]
sv.water$dap = day(sv.water$timestamp) + 4
sv.water = sv.water[sv.water$dap > 14,]
Plot the S. viridis water volume data
water.plot = ggplot(sv.water, aes(y = water.amount, x = dap,
group = factor(interaction(early.late,
treatment,dap)),
fill = factor(interaction(early.late,
treatment)),
color = factor(interaction(early.late,
treatment)))) +
geom_boxplot(outlier.colour = NULL) +
scale_x_continuous("Days after planting") +
scale_y_continuous("Water volume (ml)") +
theme_bw() +
theme(legend.position = c(0.2,0.8),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(fill='Water treatment', color='Water treatment')
water.plot = water.plot + labs(title='Figure 1B')
print(water.plot)
Calculate WUE for S. viridis and S. italica
wue.parents = function(water, biomass) {
# WUE data vectors
dap.list = c()
plant.list = c()
water.list = c()
biomass.list = c()
treatment.list = c()
group.list=c()
# Get unique barcodes for biomass
barcodes = unique(biomass[(biomass$genotype=='A10' | biomass$genotype=='B100') &
(biomass$treatment == 100 |
biomass$treatment == 33),]$plant_id)
for(barcode in barcodes) {
snapshots = biomass[biomass$plant_id==barcode,]
snapshots = snapshots[with(snapshots, order(dap)),]
for(row in 1:nrow(snapshots)) {
total.water = sum(water[water$dap <= snapshots[row,]$dap &
water$plant.barcode == barcode,]$water.amount)
dap.list = c(dap.list, snapshots[row,]$dap)
plant.list = c(plant.list, barcode)
water.list = c(water.list, total.water)
biomass.list = c(biomass.list, snapshots[row,]$fw_biomass)
treatment.list = c(treatment.list, snapshots[row,]$treatment)
group.list = c(group.list,snapshots[row,]$group)
}
}
wue.data = data.frame(plant_id=plant.list,
dap=dap.list,
water=water.list,
biomass=biomass.list,
treatment=treatment.list,
group=group.list)
return(wue.data)
}
parents.wue = wue.parents(water.data, vis.data)
parents.wue = parents.wue[parents.wue$water != 0,]
Plot S. viridis and S. italica WUE in mg/mL.
wue.plot = ggplot(parents.wue, aes(x=dap, y=(biomass/water)*1000,
color=factor(group))) +
geom_point(size=2.5) +
geom_smooth(method="loess", size=1) +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(lim=c(-2.1, 26),
name="Water-use efficiency (mg/mL)") +
theme_bw() +
theme(legend.position=c(0.2,0.8),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color='Genotype-treatment')
wue.plot = wue.plot + labs(title='Figure 4D')
print(wue.plot)
Genotype WUE function
wue = function(water, biomass, genotype) {
# WUE data vectors
dap.list = c()
plant.list = c()
water.list = c()
biomass.list = c()
treatment.list = c()
# Get unique barcodes for biomass
barcodes = unique(biomass[biomass$genotype==genotype &
(biomass$treatment == 100 |
biomass$treatment == 33),]$plant_id)
for(barcode in barcodes) {
snapshots = biomass[biomass$plant_id==barcode,]
snapshots = snapshots[with(snapshots, order(dap)),]
for(row in 1:nrow(snapshots)) {
total.water = sum(water[water$dap <= snapshots[row,]$dap &
water$plant.barcode == barcode,]$water.amount)
dap.list = c(dap.list, snapshots[row,]$dap)
plant.list = c(plant.list, barcode)
water.list = c(water.list, total.water)
biomass.list = c(biomass.list, snapshots[row,]$fw_biomass)
treatment.list = c(treatment.list, snapshots[row,]$treatment)
}
}
wue.data = data.frame(plant_id=plant.list,
dap=dap.list,
water=water.list,
biomass=biomass.list,
treatment=treatment.list)
return(wue.data)
}
Analyse pair-wise treatment differences for two-day increments for each genotype.
analyze_wue = function(wue.data) {
days = c()
diff.low = c()
diff.up = c()
pvals = c()
for(day in levels(factor(as.integer(wue.data$dap)))) {
day = as.integer(day)
control = wue.data[(as.integer(wue.data$dap) == day |
as.integer(wue.data$dap) == day + 1) &
wue.data$treatment == 100,]
drought = wue.data[(as.integer(wue.data$dap) == day |
as.integer(wue.data$dap) == day + 1) &
wue.data$treatment == 33,]
control.wue = control$biomass / control$water
drought.wue = drought$biomass / drought$water
test = t.test(x=control.wue,y=drought.wue)
days = c(days,day)
diff.low = c(diff.low, test$conf.int[1])
diff.up = c(diff.up, test$conf.int[2])
pvals = c(pvals, test$p.value)
}
results = data.frame(dap=as.numeric(days),
conf.int.low=diff.low,
conf.int.up=diff.up,
pvalue=pvals)
return(results)
}
Analyze WUE for each genotype and aggregate results
wue.results = data.frame()
for(genotype in c('A10', 'B100')) {
genotype.wue = wue(water.data, vis.data, genotype)
genotype.wue = genotype.wue[genotype.wue$water != 0,]
genotype.wue.results = analyze_wue(genotype.wue)
genotype.wue.results$genotype = genotype
wue.results = rbind(wue.results, genotype.wue.results)
}
Control for multiple testing by controlling the FDR
qvalues.wue = qvalue(wue.results$pvalue)
wue.results$qvalue = qvalues.wue$qvalues
write.table(wue.results, file='wue.stats.csv', sep=',', row.names=FALSE)
print(wue.results)
## dap conf.int.low conf.int.up pvalue genotype qvalue
## 1 11 0.0001290017 6.743361e-04 0.0055453872 A10 0.04066617
## 2 12 -0.0003304552 6.531982e-04 0.5007434856 A10 0.90056523
## 3 13 -0.0004887093 6.645100e-04 0.7558400553 A10 0.95533507
## 4 14 -0.0005521785 8.912840e-04 0.6285589728 A10 0.92188649
## 5 15 -0.0006900247 1.272981e-03 0.5413573872 A10 0.91614327
## 6 16 -0.0013750018 1.352784e-03 0.9863361510 A10 0.98633615
## 7 17 -0.0032500144 1.960050e-03 0.5946773334 A10 0.92188649
## 8 18 -0.0026416988 1.024200e-03 0.3665557006 A10 0.80642254
## 9 19 -0.0025138146 1.198806e-03 0.4687644330 A10 0.90056523
## 10 20 -0.0021111216 2.216305e-03 0.9594528855 A10 0.98633615
## 11 21 -0.0013777865 3.004615e-03 0.4446854076 A10 0.90056523
## 12 22 -0.0012554051 3.263703e-03 0.3590888551 A10 0.80642254
## 13 23 -0.0029737246 4.659146e-03 0.6177633347 A10 0.92188649
## 14 25 -0.0016604803 2.491094e-03 0.6763471555 A10 0.95533507
## 15 26 -0.0016183794 2.333780e-03 0.7068814380 A10 0.95533507
## 16 27 -0.0014579106 1.859317e-03 0.8033499452 A10 0.95533507
## 17 28 -0.0015730725 1.687446e-03 0.9422861297 A10 0.98633615
## 18 29 -0.0016425521 1.397919e-03 0.8687331415 A10 0.98633615
## 19 30 -0.0017235963 1.223038e-03 0.7270761809 A10 0.95533507
## 20 31 -0.0017633197 1.027390e-03 0.5898632801 A10 0.92188649
## 21 32 -0.0018009970 9.249981e-04 0.5116847871 A10 0.90056523
## 22 33 -0.0019234072 2.482920e-03 0.7841395855 A10 0.95533507
## 23 11 -0.0004663250 2.293599e-04 0.4748570604 B100 0.90056523
## 24 12 -0.0003827145 4.190170e-04 0.9249400490 B100 0.98633615
## 25 13 -0.0007208859 1.779908e-04 0.2235222936 B100 0.64019648
## 26 14 -0.0007089761 6.505828e-04 0.9282507680 B100 0.98633615
## 27 15 -0.0013128372 2.973716e-04 0.1995779600 B100 0.62724502
## 28 16 -0.0014265056 1.455777e-03 0.9826344845 B100 0.98633615
## 29 17 -0.0054259745 -3.416845e-04 0.0302482577 B100 0.12099303
## 30 18 -0.0039584923 7.814527e-05 0.0582767730 B100 0.21368150
## 31 19 -0.0053958070 -8.822144e-04 0.0097441669 B100 0.06124905
## 32 20 -0.0049249253 -5.837513e-04 0.0161879273 B100 0.07914098
## 33 21 -0.0062171326 -1.283934e-03 0.0055344094 B100 0.04066617
## 34 22 -0.0056491252 -1.224773e-03 0.0043669943 B100 0.04066617
## 35 23 -0.0061377032 -1.581340e-03 0.0036480639 B100 0.04066617
## 36 25 -0.0056519496 -1.646284e-03 0.0009985777 B100 0.03550789
## 37 26 -0.0050694636 -1.330337e-03 0.0016139948 B100 0.03550789
## 38 27 -0.0041897428 -5.455619e-04 0.0130667633 B100 0.07186720
## 39 28 -0.0036907493 -2.531025e-04 0.0261291938 B100 0.11496845
## 40 29 -0.0033005146 1.619285e-04 0.0734489761 B100 0.24859653
## 41 30 -0.0026887815 6.877754e-04 0.2327987196 B100 0.64019648
## 42 31 -0.0025257356 7.093671e-04 0.2546041527 B100 0.65897545
## 43 32 -0.0024342851 8.292998e-04 0.3164971519 B100 0.77365970
## 44 33 -0.0028628094 2.203841e-03 0.7831321502 B100 0.95533507
In this section we use linear modeling to estimate tiller number.
Calculate height-width ratio.
vis.data$height_width_ratio = vis.data$height_above_bound / vis.data$extent_x
Download data for manually measured tiller counts 195 random images.
if (!file.exists('200tiller_counts.txt')) {
download.file('http://files.figshare.com/1845359/200tiller_counts.txt',
'200tiller_counts.txt')
}
Read manual tiller count data.
tiller.counts = read.table(file="200tiller_counts.txt",sep="\t",header=TRUE)
Format date.
tiller.counts$date = as.POSIXct(paste(paste(tiller.counts$year,
tiller.counts$month,
tiller.counts$day, sep='-'),
paste(tiller.counts$hours,
tiller.counts$minutes,
tiller.counts$seconds, sep=':'), sep=' '))
Merge with VIS snapshot table.
tiller.counts = merge(x=tiller.counts, y=vis.data, by = 'date')
Model tiller counts.
tiller.model.full = lm(ave_tillers ~ fw_biomass + height_above_bound + solidity +
extent_x + height_width_ratio, data=tiller.counts)
Variance inflation factor.
vif(tiller.model.full)
## fw_biomass height_above_bound solidity
## 17.338117 22.493335 1.651753
## extent_x height_width_ratio
## 14.767073 3.280364
Drop height and recalculate VIF
tiller.model1 = lm(ave_tillers ~ fw_biomass + solidity + extent_x +
height_width_ratio, data=tiller.counts)
vif(tiller.model1)
## fw_biomass solidity extent_x
## 10.145172 1.393258 11.362312
## height_width_ratio
## 1.855525
Drop extent x.
tiller.model2 = lm(ave_tillers ~ fw_biomass + solidity + height_width_ratio,
data=tiller.counts)
vif(tiller.model2)
## fw_biomass solidity height_width_ratio
## 1.027223 1.044973 1.039641
summary(tiller.model2)
##
## Call:
## lm(formula = ave_tillers ~ fw_biomass + solidity + height_width_ratio,
## data = tiller.counts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.7006 -1.2510 -0.0308 0.8410 5.6311
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.17731 0.54049 9.579 < 2e-16 ***
## fw_biomass 0.22888 0.01275 17.957 < 2e-16 ***
## solidity -0.07202 1.87918 -0.038 0.969
## height_width_ratio -2.16290 0.32335 -6.689 2.42e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.699 on 191 degrees of freedom
## Multiple R-squared: 0.6501, Adjusted R-squared: 0.6446
## F-statistic: 118.3 on 3 and 191 DF, p-value: < 2.2e-16
Drop solidity
tiller.model3 = lm(ave_tillers ~ fw_biomass + height_width_ratio, data=tiller.counts)
summary(tiller.model3)
##
## Call:
## lm(formula = ave_tillers ~ fw_biomass + height_width_ratio, data = tiller.counts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.7003 -1.2525 -0.0264 0.8392 5.6326
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.1640 0.4136 12.487 < 2e-16 ***
## fw_biomass 0.2289 0.0126 18.173 < 2e-16 ***
## height_width_ratio -2.1650 0.3177 -6.815 1.18e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.694 on 192 degrees of freedom
## Multiple R-squared: 0.6501, Adjusted R-squared: 0.6465
## F-statistic: 178.4 on 2 and 192 DF, p-value: < 2.2e-16
Plot partial regressions.
avPlots(model = tiller.model3, pch=16, layout=c(2,1), main="Figure 5B")
Download data for manually measured tiller counts for 646 random images.
if (!file.exists('660tiller_counts.txt')) {
download.file('http://files.figshare.com/1845358/660tiller_counts.txt',
'660tiller_counts.txt')
}
Predict tiller count for a second set of ~660 randomly selected plants.
tiller660 = read.table(file="660tiller_counts.txt", sep="\t", header=TRUE)
tiller660 = merge(x=tiller660, y=vis.data, by='datetime')
tiller660$predicted_tiller_count = predict.lm(object = tiller.model3,
newdata=tiller660)
Plot manual versus predicted tiller counts.
tiller.plot = ggplot(data=tiller660, aes(x=tiller_count, y=predicted_tiller_count,
color=factor(genotype))) +
geom_point(size=2.5) +
geom_abline(intercept=0, slope=1) +
scale_x_continuous(lim=c(-1,14), "Manual tiller count") +
scale_y_continuous(lim=c(-1,14), "Predicted tiller count") +
theme_bw() +
theme(legend.position=c(0.45,0.95),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold"),
legend.direction="horizontal") +
labs(color="Genotype")
tiller.plot = tiller.plot + labs(title='Figure 5C')
print(tiller.plot)
Predict tiller counts for the whole data set.
vis.data$tiller_count = predict.lm(object = tiller.model3, newdata = vis.data)
Plot tiller counts for S. viridis and S. italica water treatments 100% and 33% full-capacity.
tc.par.plot = ggplot(vis.data[(vis.data$genotype == 'A10' |
vis.data$genotype == 'B100') &
(vis.data$treatment == 100 |
vis.data$treatment == 33),],
aes(x=dap, y=tiller_count, color=factor(group))) +
geom_point(size=2.5) +
geom_smooth(method="loess",size=1) +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(lim=c(0,12),
name="Estimated tiller count") +
theme_bw() +
theme(legend.position=c(0.25,0.75),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Genotype-treatment")
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 2 rows containing missing values (geom_point).
tc.par.plot = tc.par.plot + labs(title='Figure 5D')
print(tc.par.plot)
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 2 rows containing missing values (geom_point).
Plot estimated tiller counts for all genotypes.
tiller.facet.plot = ggplot(vis.data[vis.data$treatment == 100 |
vis.data$treatment == 33,],
aes(x=dap, y=tiller_count, color=factor(treatment))) +
geom_point(size=1) +
geom_smooth(method="loess",size=1) +
scale_x_continuous(name="Days after planting") +
scale_y_continuous(lim=c(0,12),
name="Estimated tiller count") +
theme_bw() +
theme(legend.position=c(0.8,0.1),
axis.title.x=element_text(face="bold"),
axis.title.y=element_text(face="bold")) +
labs(color="Treatment") +
facet_wrap(~ genotype, ncol = 3)
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 3 rows containing missing values (stat_smooth).
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 2 rows containing missing values (geom_point).
## Warning: Removed 2 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 4 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_point).
tiller.facet.plot = tiller.facet.plot + labs(title='Supplemental Figure S4')
print(tiller.facet.plot)
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 3 rows containing missing values (stat_smooth).
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 1 rows containing missing values (stat_smooth).
## Warning: Removed 2 rows containing missing values (geom_point).
## Warning: Removed 2 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 4 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_point).
Output a table of all of the traits input and added here. Export all of the VIS traits (except color).
h.table = vis.data
h.table$plant_id = as.character(h.table$plant_id)
h.table$wue = NA
for(r in 1:nrow(h.table)) {
row = h.table[r,]
total.water = sum(water.data[water.data$dap <= row$dap & water.data$plant.barcode == row$plant_id,]$water.amount)
h.table[r,]$wue = row$sv_area / total.water
}
write.table(h.table,file="vis.traits.csv",quote=FALSE,sep=",",row.names=FALSE)